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Abstract 

We report spontaneous translocation of small interfering RNA (siRNA) inside carbon nanotubes (CNTs) 
of various diameters and chirality using all atom molecular dynamics (MD) simulations with explicit sol- 
vent. We use Umbrella sampling method to calculate the free energy landscape of the siRNA entry and 
translocation event. Free energy profiles shows that siRNA gains free energy while translocating inside 
CNT and barrier for siRNA exit from CNT ranges from 40 to 1 10 kcal/mol depending on CNT chirality 
and salt concentration. The translocation time x decreases with the increase of CNT diameter with a critical 
diameter of 24 A for the translocation. In contrast, double strand DNA (dsDNA) of the same sequence 
does not translocate inside CNT due to large free energy barrier for the translocation. This study helps in 
understanding the nucleic acid transport through nanopores at microscopic level and may help designing 
carbon nanotube based sensor for siRNA. 
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I. INTRODUCTION 



Small interfering RNA (siRNA) is 21-23 nucleotides long double stranded RNA (dsRNA) 
which has been demonstrated to treat diseases like cancer and hepatitis. For such medical 
applications, siRNA has to be delivered inside the cell without degradation by complexing it 
with delivering vectors. However, efficient method for the delivery of siRNA to the target cell is 
lacking. Carbon nanotubes (CNTs) are cylindrical rolled graphene sheets which have emerged 
as new nano material to be used for biomedical applications due to their potential strengths 
[1-4]. CNTs have ultrahigh surface area which make them suitable drug delivery vectors through 
specific binding of biomolecules, drugs on their outer surface. This unique feature along with the 
other physical, thermal and optical properties enable high throughput drug delivery and medical 
imaging [2]. CNTs are not soluble in many solvents and also known to be toxic to the cell. 
This problem has been addressed by covalent or non-covalent surface functionalization of CNTs 
[1-4]. Delivery of RNA by translocating RNA-CNT hybrid through membranes of MCF7 breast 
cancer cells has also been studied [5]. siRNA delivery through pristine and functionalized CNT 
have been investigated by several groups in recent years. In this delivery approach siRNA is 
adsorbed on the CNT. After transfecting siRNA-CNT hybrid inside cells, the siRNA dissociates 
from the CNT. Recently we have shown that siRNA can adsorb onto the surface of the CNT and 
graphene by unzipping its base-pairs and binds strongly with high stability which can be used to 
deliver siRNA [6, 7]. All these studies are concerned with siRNA adsorption/binding on CNT 
surface only. In this paper, we propose encapsulation of siRNA inside CNT as a possible delivery 
mechanisms. 

Transport of biomolecules across cell membranes and nanopores is an important process in 
living organisms [8]. Single stranded DNA (ssDNA) transport in biologically occurred nanometer 
scale cc-hemolysin ion channel has been observed by Kasianowicz et. al. [9] in 1996. They 
observed that the translocation of ssDNA depends on the DNA concentration, applied voltage 
and structure of the nanopore. Since then, biomolecular transport in nanopores has become a 
major research activity both from fundamental and application point of view [10-12]. Polymer 
translocation through membranes [13] and nanopores [11, 14, 15] has been studied theoretically 
as a diffusion process across free energy barrier that arises due the chemical potential gradients. 
It has been found that the translocation time is proportional to the polymer length. Sophisticated 
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nanopore based DNA sequencing techniques are now-a-days feasible to determine the sequence of 
DNA by measuring ionic currents as ssDNA translocates through the nanopore [16, 17]. However, 
microscopic level details are necessary for the better understanding of the DNA sequencing at the 
single base-pair level which is the ultimate goal of DNA sequencing technology. Translocation 
of ions, nucleic acids through solid state nanopores have been studied by Dekker and co-workers 
[16, 18, 19]. Translocation of oligonucleotides inside CNT have been reported as well [20-25]. 
Recently, translocation of DNA through nanopores in suspended graphene sheets have also been 
studied experimentally that finds applications in DNA sequencing [26-29]. 

Based on our MD simulation studies, we propose two different methods of siRNA delivery 
using CNTs and graphene. The first mechanism is based on the strong adsorption of siRNA on the 
CNT/graphene surface [6, 7] and the second one is by encapsulation of siRNA through translo- 
cation inside the CNT. The translocation of siRNA inside a CNT has not been studied earlier. 
In this paper, we report for the first time, the translocation studies of siRNA inside CNT. This 
translocation is driven by the favorable van der Waals (vdW) interaction between the siRNA and 
CNT. To understand the free energy landscape of siRNA translocation inside CNT, we have also 
calculated the free energy profile using Umbrella sampling method [30-32] for CNTs of different 
chirality and at different salt concentrations. The free energy profile indicates that siRNA gains en- 
ergy inside CNT and faces a large energy barrier to escape from the interior of CNT. Insights into 
siRNA-CNT interaction and structural changes will help in understanding the microscopic picture 
of the siRNA delivery mechanism by CNTs. The translocation time x of siRNA is decreasing with 
increasing the diameter of the CNT. Our simulation results can guide experimental approach for 
siRNA translocation inside CNT for sequencing and drug delivery applications as well. The rest 
of the of the paper is organized as follows: in the section II, we give the details of the simulation 
methods, followed by the discussion on the results of spontaneous translocation. In section III, 
we discuss the results on free energy landscapes using MM-GBSA as well as Umbrella sampling 
method. Finally in section IV, we give a summary of the main results and conclude. 

II. METHODS 

We have built single walled CNT (SWCNT) of different chiralities and diameters of length 
approximately 87 A. The diameters d of the built SWCNTs are 27.01 A, 26.74 A, 25.40 A, 24.06 
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A and 22.73 A for (25, 15), (20, 20), (19, 19), (18, 18) and (17, 17) nanotubes respectively. We 
study several systems for nucleic acid translocation inside CNTs: (i) siRNA translocation mech- 
anism inside (20, 20) CNT at and 107 mM NaCl concentration to investigate the effect of salt 
concentration (ii) siRNA translocation inside (19, 19), (18, 18) and (17, 17) CNTs to study the 
effect of diameter (iii) siRNA translocation mechanism inside (25, 15) to study the effect of chi- 
rality and finally (v) dsDNA (of identical sequence except adenine instead of uracil of siRNA) 
translocation mechanism inside (20, 20) CNT to understand the difference between the interaction 
of siRNA and dsDNA with CNT. For all these cases, the siRNA/dsDNA was placed close to one 
of the open ends of CNT such that the CNT axis and the siRNA axis are parallel. We have used 
the same siRNA sequence and structure as described in our earlier publication [6, 7]. We have 
used ff99 force field [33] for all the simulations reported in this paper. The siRNA-CNT system 
is solvated with water and counterions. To make a charge neutral system for different salt concen- 
tration, appropriate number of Na + and Cl~ ions were added. The system was then subjected to 
standard simulation protocol [34, 35] with periodic boundary conditions applied in all three direc- 
tions. Non-bonded interactions are truncated at 9 A and the long range electrostatic interactions 
were calculated with the Particle Mesh Ewald (PME) method [36]. Finally for trajectory analysis, 
simulations are performed in constant volume-constant temperature (NVT) ensemble up to 100 
ns depending on the diameter of CNT. For the Umbrella sampling method, 120-130 ns of NVT 
simulations are performed for each system. Full system details such as number of NaCl ions, wa- 
ter residues, box dimensions and total number of atoms are summarized in Table I. To check the 
effect of force field, we have also performed the simulation with (20, 20) CNT and siRNA using 
ff 10 force field [37]. 

IH. RESULTS AND DISCUSSION 

A. Translocation of siRNA inside CNTs 

Snapshots shown in Figure 1 represent the translocation process of siRNA inside (20, 20) CNT 
at various instants of time. To present the details of the system studied, we show the initial siRNA- 
(20, 20) CNT system with counterions and water in Figure 1(a). In Figures 1(b) to l(i), we show 
the instantaneous snapshots of the system in few ns interval. For clarity, water and counterions are 
not shown. From the snapshots and subsequent analysis, we observe that the siRNA translocates 



4 



into the interior of (20, 20) CNT from the initial configuration. The translocation of siRNA inside 
(20, 20) CNT is very fast and happens in 17.4 ns. After the translocation, siRNA stays inside 
CNT with least positional fluctuations during the remaining simulation period of 30-40 ns. We 
have also observed the translocation of siRNA at 107 mM, inside (19, 19), (18, 18) and (25, 15) 
CNTs. The siRNA inside (20, 20) CNT at mM has least structural deformations compared 
to all other CNTs of different diameter and chirality studied in this work. Similar results are 
observed at 107 mM NaCl concentration also. In all these cases, the siRNA is very stable inside 
CNT without any further movement after the translocation. This siRNA-CNT hybrid can serve 
as siRNA delivery vector for the delivery of siRNA for RNAi therapy. For the siRNA release 
after delivery, functionalizing siRNA with nanoparticles or application of electric field can be used. 

To understand the energetics of siRNA as it translocates inside (20, 20) CNT, we have calcu- 
lated the van der Waals (vdW) interaction energy <|>(r) between siRNA and CNT at various instants 
of time. We plot <|>(r) as a function of the siRNA spatial position r at various instants of time t in 
Figure 2. Figure 2 shows how different parts of siRNA interact with CNT and also the interaction 
strength along CNT axis h. In Figure 2, r = corresponds to the center of mass (COM) position 
of siRNA. For the calculation of <|>(r), siRNA has been divided into 4 parts along its helix axis. 
The data shown in Figure 2 represents how the COM of these 4 parts of siRNA interact with CNT 
as it translocates. It can be seen from the Figure 2 that <|>(r) has different shapes at t = ns, 10 
ns, 15 ns, 25 ns and 40 ns. Initially siRNA is outside CNT with only two sticky-ends (these are 
un-paired nucleobases in a nucleic acid) and one intact Watson-Crick (WC) base-pair lying inside 
CNT. Hence the interaction of siRNA with CNT is very less which are far apart from each other 
(see t = ns curve in Figure 2). As t increases, these sticky-ends interact strongly with CNT as 
shown in Figure 2 at t = 10 ns, 15 ns, 25 ns and 40 ns. The siRNA strongly interacts with CNT 
after complete translocation for t > x; where x is the translocation time. The increase of § after 
complete translocation is about 407 kcal/mol or 690.8 ksT. 



5 



B. Free energy landscape during translocation 



1. Free energy using equilibrium simulation: MM-GBSA and 2PT method 

In general, the binding free energy for the non-covalent association of two molecules A and 
B may be written as AG (A + B—> AB) — Gab — Ga — Gb- For any species on the right hand side 
G(X) = H(X) — TS(X), where H is the enthalpy, T is the absolute temperature and S is the entropy 
of the molecule. Therefore the binding energy at constant temperature can be written as 

AG = AH -T AS (1) 

The calculation of enthalpy difference (AH) and entropy difference (AS) contributions to the 
binding free energy in Eqn. 1 are done using molecular mechanics - generalized Born surface 
area (MM-GBSA) method and 2 phase thermodynamic (2PT) method [38-40], respectively. The 
change in enthalpy AH can be decomposed into the gas-phase energy (AE gas ) and solvation free 
energy (AG so i); i.e., AH = AE gas + AG so i. For any species A, B or complex AB, the gas-phase 
energy, E gRS is calculated from molecular mechanics using E gas = E e \ e + + E- mt , where, E e \ e 
is the electrostatic energy calculated from the Coulomb potential, £Vd w is the non-bonded van der 
Waals energy and E{ nt is the internal energy contribution arising from bond stretching, angle bend- 
ing and dihedrals. On the other hand, the solvation free energy, G so i can be further decomposed as 
electrostatic G e i and non-electrostatic G non -ei contributions; i.e., G so i = G e i + G non -ei- The electro- 
static energy, G e i is calculated from Generalized Born (GB) method which assumes that the atoms 
in a molecule are spheres of radius 7?, (called Born radius) and have an effective charge q t (called 
Born charge). The molecule is assumed to be surrounded by a solvent of dielectric constant £ (80 
for water at 300 K) and the solute atoms have a dielectric constant of 1. The analytic expression 
for the G e i in the GB model [41, 42] is given by 

1 y m± 

61 " ijjfGBinj^Rj) 

where 

f GB {r ih R u Rj) = [rfj +R l R j cxp {-r^/AR^)] * (3) 

where nj is the distance between atoms /, j and i?,- are effective Born radii and K is the Debye- 
Hiickel screening parameter. The non-electrostatic energy, G non -ei is calculated as yx SASA + (3; 



1 exp(-K/ GB ) \ 



(2) 
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where y is the surface tension parameter (y = 0.0072 kcal/mol-A 2 ; (3 = kcal/mol) and SASA is 
the solvent-accessible surface area of the molecule. The entropy contribution appearing in Eqn. 1 
is carried out using 2 phase thermodynamic approach proposed by Lin et. al. [38-40] which is 
motivated by the observation that the density of states (S(v)) of a liquid can be decomposed into 
a gas component and a solid component [38-40]. The density of states function can be calculated 
from the Fourier transform of velocity auto-correlation function; i.e., 

S(y) = 7^- lim T V V m;cHt)e- l2%vt dt (4) 

where cj k (t) is the velocity auto-correlation function which is given by 

Cj\t) = lim 1 f Vj k (t) Vj k (t + t')dt'. (5) 

ZT J—z 

The entropy S is then calculated from the knowledge of S(v) using 

S = k B \nQ + $- ld ^ (6) 

Ol N,V 

where the partition function Q is given by 

poo 

lnQ= / dvS(v)lnq H0 (v). (7) 
Jo 

The decomposition of 5(v) into gas-phase and solid-phase is performed to get the entropy of 
molecular species A, B or complex AB to get finally the change in free energy TAS. The 2PT 
method has found successful application in several related problems [38-40, 43-45]. 



The binding free energy (AG) as function of time is shown in Figure 3(a). Entropy (TAS) 
and enthalpy (AH) contributions to AG are shown in the inset. AS decreases with time since the 
fluctuations in siRNA nucleobases are suppressed by CNT making less microstates available for 
siRNA as translocation progresses. However, AH is increasing with time and attains constant 
value after complete translocation. AG also follows similar trend since the entropy contribution 
is very small compared to enthalpy contribution. We plot the vdW contribution A§ to the total 
binding energy in Figure 3(b). When siRNA and CNT are far apart in the initial stage, A(|) is very 
less but increases as siRNA translocates inside CNT. We also show instantaneous snapshots of the 
siRNA-CNT hybrid at various instants of time t in Figure 3(b). A<|) is -900 kcal/mol after complete 
translocation of siRNA compared to -85 kcal/mol at the initial stage. The interaction between 
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siRNA and CNT is driven mainly by vdW interaction (Figure 3(b)). This is similar to our earlier 
results that the unzipping and adsorption of siRNA on CNT/graphene is also driven by the vdW 
interaction [6, 7]. The conversion of time into distance between CNT and siRNA while plotting 
AG(r) has been discussed in the supplementary materials [46]. 

2. Free energy using Umbrella sampling method 

In the un-biased MD simulations discussed so far, siRNA has translational motion in 1-D along 
nanotube axis h in order to translocate into the interior of CNT. The free energy landscape of the 
siRNA translocation process along such suitable reaction coordinate will help in understanding 
the stable and unstable states for siRNA with respect to the CNT. In the previous section we 
have described the free energy calculation from the equilibrium un-biased MD simulation using a 
combination of MM-GBSA and 2PT methods. Such methods in equilibrium MD cannot properly 
account for the sampling of high energy states. For this reason, we use Umbrella sampling 
method [30] to sample the entire phase space along chosen reaction coordinate and calculate 
the free energy profile for siRNA translocation inside CNT. The basic idea implemented in 
constructing the free energy is to add a biasing potential to obtain sampling of less probable states 
for siRNA translocation process and record the biased histograms. This is done for a series of 
biasing potentials at various values of chosen reaction coordinate that span the entire translocation 
path of interest. We make sure that the successive histograms have enough overlap in order to 
reconstruct the un-biased free energy from biased histograms. Below we give a brief overview of 
the Umbrella sampling method. For further details, readers are refereed to the excellent text book 
by Frenkel and Smit [32]. 

With the unperturbed potential V(r^), we have added a biasing harmonic potential t/;(£) = 
Ikfe-^i) 2 , resulting in the perturbed potential V'^) = V(r N ) + 1 £f= 1 */f(£); where is the re- 
strained distance between the center of masses of 'far end' of CNT and first two Watson-Crick 
hydrogen bonded base-pairs of siRNA close to the CNT for z' th window. The values of k and A£, 
should be chosen optimally such that the phase space is properly sampled within reasonable time 
scales. Therefore, we have optimized the force constant k of A£, and equilibration time 

before performing the Umbrella sampling. The optimized values are k = 4 kcal/mol-A 2 , A^ ; = 1 
A and equilibration time of 1 ns to ensure proper sampling and overlapping of successive his- 
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tograms, The values of ^ is changed from 130 A to 10 A until the complete translocation 

happens in steps of 1 A totaling N W = 120 Umbrella simulations (windows) of each 1 nano second 
duration. The probability distribution of the reaction coordinate £ separating CNT and siRNA is, 

Similarly the probability distribution in the presence of is 

P,%) = j, J dr»S&-t,)e-^^>\ 

Re-arranging this gives, 

m) = ^ m) Pi'(t>) (8) 

here, 8 is a Dirac delta function, N is the number of atoms, denotes the set of atom coordinates, 
(3=1 /kgT, where kg is the Boltzmann constant, T is the absolute temperature, V(rj is the po- 
tential energy and Z, Z' are the partition functions of un-biased and biased systems, respectively. 
Here problem arises due to the fact that P(£) becomes exceedingly small for values of £ which give 
significant contribution to the free energy. Umbrella sampling makes use of a biasing potential to 
sample the region of phase space for which P(£) is exceedingly small. Thus, the distribution func- 
tion P(£) can be obtained using Eqn. 8 (to within a multiplicative constant) from the measurement 
of the biased distribution P'(Q- The un-biased free energy (F) was constructed self-consistently 
using the weighted histogram analysis method (WHAM) [30-32] with the following Eqns. 9 and 
10. 

p(n = jw=Wj(j) (q\ 

where F\ is given by 

e-M = J d&-W®P{%) (10) 

where n,- is the number of data points in z' th window. By piecing together the relative free energies 
measured using a number of biasing potentials, we construct F(Q over the chosen range of £, = 
10 A to 130 A. 



Figure 4(a) shows the free energy for siRNA in (20, 20) CNT at mM NaCl concentra- 
tion. We find that when siRNA is outside (20, 20) CNT, the free energy is zero since there is no 
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interaction between siRNA and (20, 20) CNT. The snapshots at various time instants show the 
translocation process as £ is decreasing from 130 A to 10 A. For £ > 105 A, siRNA is outside 
(20, 20) CNT and does not interact with it. As the value of £, is decreased from 105 A, F(^) 
decreases to minima at 92 A at which siRNA is in a local stable state. At this stage siRNA rotates 
inside CNT and two unpaired bases of siRNA and one WC base-pair unzips to interact with (20, 
20) CNT as can be seen in Figure 4(a). Further decrease in £, results in increasing F(^) having 
an energy barrier for siRNA to translocate inside CNT. F(^) reaches maximum at ^ = 80 A with 
an energy barrier of AG = 9 kcal/mol. The energy barrier for further translocation arises due to 
strong vdW interaction between unpaired nucleobases of siRNA and CNT. Thermal fluctuations 
help siRNA to overcome this barrier. Note that, once it overcomes this barrier, favorable vdW 
interaction between nucleobases and CNT helps siRNA to translocate inside. From ^ = 80 A, 
siRNA goes inside (20, 20) CNT easily and F(^) decreases until complete siRNA translocation at 
^min = 42 A. This is the most stable position for siRNA inside (20, 20) CNT where siRNA stays 
for the rest of the simulation time. 

Figure 4(b) presents free energy for siRNA in CNTs for (20, 20), (25, 15) chiralities at 
mM NaCl concentration, for (20, 20) CNT at 107 mM NaCl concentration and for dsDNA in 
(20, 20) CNT at mM NaCl concentration. In the un-biased simulation, we have observed the 
translocation of siRNA in CNTs but no translocation of dsDNA is observed. The dsDNA has 
favorable state only outside the (20, 20) CNT as can be seen in Figure 4(b). This might be due 
to the less favorable interaction of thymidine nucleobase with nanotube compared to the uridine 
nucleobase interaction with nanotube [6, 7]. For ^ < L/2, where L is the length of the CNT, F(^) 
is constantly increasing which means dsDNA encounters unfavorable free energy landscape when 
translocated inside (20, 20) CNT. This dramatic difference of interaction between CNT and siRNA 
and CNT and dsDNA may be due to the relatively weaker interaction strength of thymidine with 
CNT than that of uridine with CNT [6, 7]. Moreover, the stronger WC base-pairing interaction 
energy of A-T compared to A-U [6, 7, 47, 48] makes dsDNA difficult to get unzipped. Earlier we 
have shown that unzipped base-pairs facilitate binding with CNT/graphene [6, 7]. dsDNA requires 
more than 100 kcal/mol energy to overcome a free energy barrier to translocate inside (20, 20) 
CNT which is not possible without any external force. However, possibility of the translocation 
of dsDNA in CNT of large diameter has not been studied in our simulation. In experiments of 
dsDNA translocation inside CNT, the diameter of CNT is 50-100 nm [21, 49]. In the experiments 
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of ssDNA translocation through a- hemolysin ion channel [9], the diameter of a- hemolysin is 26 
A and the polymer is single stranded. Diameter of (20, 20) CNT is 26.74 A which is close to 
the diameter of oc-hemolysin ion channel and so the channel diameter is not sufficient for dsDNA 
translocation. We have observed the translocation of ssDNA through CNT of diameter 26.74 
A [50]. Our results on ssDNA and dsDNA translocation are consistent with experiments [9]. 
F(Q for dsDNA has global minimum only outside the (20, 20) CNT. When dsDNA is forced to 
translocate inside (20, 20) CNT with an external Umbrella potential, dsDNA deforms largely by 
breaking most of the WC base-pairs due to the large external force by Umbrella potential. To have 
further confirmation that siRNA translocate inside CNT and dsDNA does not, we have performed 
two separate simulations of siRNA and dsDNA by keeping them initially inside the (20, 20) CNT. 
As expected, siRNA stays inside the CNT for long time where as dsDNA comes out of the CNT 
within 500 ps. These results confirm that siRNA spontaneously translocates inside a (20, 20) CNT 
without any external force where as dsDNA requires an external force in order to be translocated 
inside (20, 20) CNT. 

To test the effect of chirality as well as salt concentration on the translocation event, we have 
calculated the free energy profile of siRNA translocation inside (25, 15) CNT at mM NaCl 
concentration and for (20, 20) CNT at 107 mM NaCl concentration using Umbrella sampling 
method. The free energy profile for these cases are shown in Figure 4(b). In both these cases, the 
free energy minima is lower compared to the case of siRNA in (20, 20) CNT at mM concentration 
of NaCl. So for these cases we expect translocation to happen faster compared to the case for 
siRNA in (20, 20) CNT at mM NaCl concentration. Translocation time for various cases has 
been discussed in section IIIC. The inset of Figure 4(b) shows the minima of F(^). We observe 
that ^ m i n = 37 A at 107 mM as compared to ^ m i n = 42 A at mM. However, ^ m i n = 21 A for 
siRNA in (25, 15) CNT. The origin for this observation is discussed in the next section by analyzing 
un-biased simulations. 

C. Translocation time (x) 

The translocation of siRNA strongly depends on the diameter d of CNT since the vdW interac- 
tion between siRNA and CNT increases as d is decreased (nearly comparable to siRNA diameter). 
The translocation time, x is defined as the time required for at least half of the siRNA base-pairs 
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to enter inside the CNT; i.e., r CO m 5; L/2, which corresponds to f com < 43.5 A or £ < 59 A, where 
L = 87 A is the length of the CNT. We plot 

/"com and x in Figures 5(a) and 5(b), respectively 
for various diameters of CNT. In Figure 5(a), we see that r com is decreasing with time as siRNA 
moves towards the interior of the CNT. However, the sticky-ends and unzipped nucleobases of 
siRNA interact with CNT in the course of translocation, and these give rise to very rugged free 
energy profile of siRNA during the translocation. We have calculated x by tracking r com as a func- 
tion of time as shown in Figure 5(a). The dashed horizontal curve drawn at L/2 serves as the 
reference for translocation criteria. We find that x is linearly decreasing with CNT diameter d. 
Similar dependence was observed for DNA transport through graphene nanopore [26] where the 
ionic conductance was found to be proportional to the pore diameter. We note that our results on 
x are very sensitive to the initial model building. To see the effect of initial relative positioning of 
siRNA with respect to CNT, we simulated the translocation process for three different initial con- 
ditions for (20, 20) CNT. We find x to be 18 ns, 14 ns and 17.4 ns and see a strong dependence of x 
on the initial conditions. Hence our results on x can only serve as qualitative understanding of the 
systems under investigation in this study. To understand the origin of this diameter dependence, 
we have calculated the number of contacts N c between the siRNA and CNTs and have shown them 
in Figure 5(c). N c of siRNA is calculated within 5 A from the inner surface of CNT that mostly 
represents the effective vdW interaction range. As shown in Figure 5(c), N c of siRNA are 1 10 at t 
= ns and rapidly increases to a constant value after complete translocation. The maximum value 
of N c increases with d ranging from 520 inside the (18, 18) CNT to 815 inside the (25, 15) CNT. 
siRNA translocation happens above a critical CNT diameter of 24.0 A which corresponds to (18, 
18) CNT, below which no translocation is observed. For siRNA in (17, 17) CNT and dsDNA in 
(20, 20) CNT, the value of N c is constant and very less indicating no translocation in the simulated 
time scale. Since we have not observed the translocation of siRNA in (17, 17) CNT and dsDNA in 
(20, 20) CNT, x is assumed to be very large. For the translocation of siRNA inside CNT, the min- 
imum diameter of CNT, d m i n should be greater than or equal to the diameter of siRNA + effective 
vdW radius of CNT and siRNA; i.e., d m m ^ ^siRNA ^vdw • The diameter of (17, 17) CNT is 22.73 
A which is very close to siRNA diameter and hence cannot accommodate siRNA unless siRNA is 
severely stretched. As expected, with decreasing CNT diameter, the deformation in siRNA is more 
due to strong vdW interaction between siRNA and CNT. The structural aspects of siRNA during 
translocation are discussed in section HID. Interestingly, when the salt concentration is increased 
to 107 mM, x decreases drastically to 4 ns compared to 17.4 ns at mM salt concentration. High 
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salt provides better counterion condensation around phosphate atoms in siRNA backbone com- 
pared to mM case. The electrostatic screening increases as salt is increased from to 107 mM 
which reduces phosphate -phosphate electrostatic repulsion in the backbone of siRNA. Therefore, 
the stability to siRNA increases with increasing salt concentration which reduces the propensity 
of siRNA unzipping [6]. As a result, at 107 mM NaCl concentration, there is less binding of the 
siRNA with the CNT So the translocation is faster at 107 mM salt concentration and hence x is 
small. At 107 mM of NaCl concentration, the maximum value of N c is 730 which is 20 more in 
number than that of at mM. 

D. Structural aspects of siRNA during translocation inside CNT 

In order to translocate, siRNA has to cross the free energy barrier arising due to loss of 
entropy of siRNA. In this process, many structural deformations occurs in siRNA. We have 
calculated several quantities that quantify the degree of structural deformations in siRNA. In 
Figure 6(a), we we plot the root mean square deviation (RMSd) as a function of time as siRNA 
translocates in the interior of CNT for various CNT diameters as well as for different salt 
concentrations. The RMSd was calculated with respect to the initial minimized structure of 
siRNA. The RMSd of the siRNA which is inside of a thinner CNT is larger than that of the 
siRNA which is inside a fatter CNT. This difference is due to the large structural deformation 
of siRNA inside CNT of smaller diameter. The average RMSd of siRNA is 3.0 A, 5.4 A and 
8.2 A inside (20, 20), (19, 19) and (18, 18) CNT, respectively. Also RMSd for siRNA in (20, 
20) CNT is higher at 107 mM salt concentration compared to mM concentration. Among 
all the cases, siRNA in (20, 20) CNT at mM has least RMSd after translocation with small 
fluctuations. In the delivery application, the translocation method is more appropriate and di- 
ameter of the CNT can be appropriately chosen to have minimum deformation of siRNA structure. 

To translocate inside smaller diameter CNTs, siRNA has to get unzipped. To quantify the 
amount of unzipping of siRNA inside CNT, we have calculated the number of WC H-bonds 
based on geometry criteria. WC H-bonds play a crucial role in maintaining the double stranded 
form of nucleic acids and their function [51-53]. For WC H-bond calculation, we have used 
geometry based criteria wherein for H-bond, donor- acceptor (D-A) distance < 2.7 A and angle 
ZDHA > 130°. We have plotted the WC H-bonds of siRNA in various CNTs in Figure 6(b). 
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There are only three broken H-bonds out of total 48 possible H-bonds for siRNA in the case of 
(20, 20) CNT at and 107 mM salt concentrations and dsDNA in the case of (20, 20) CNT at 
mM salt concentration. As d is decreasing, more H-bonds are broken in siRNA due to stronger 
vdW interaction between siRNA and CNT. Interestingly, large number of H-bonds of siRNA are 
broken inside (25, 15) CNT. The aromatic rings of (25, 15) CNT are oriented differently with 
respect to CNT axis h. Hence nucleobases try to orient along the aromatic ring direction due 
to vdW interaction resulting in large number of broken H-bonds. The larger number of broken 
H-bonds for smaller diameter is also reflected in the larger RMSd of siRNA in case of smaller 
diameter as shown in Figure 6(a). 

Another very interesting observation is the rotation of siRNA to optimize the nucleobase ori- 
entation with respect to the aromatic rings of CNT while translocating inside CNT. We have cal- 
culated the distribution of rotation angle, P(0) of siRNA where 9 is the average rotation angle 
of siRNA as it translocates inside nanotube with respect to the structure outside nanotube. The 
rotation angle is calculated using 



where p ; (Y) is the vector joining phosphate-phosphate atoms of i th base-pair for i = 4, 5 . . . 17 at 
time t. The angular brackets (. . . ) ■ denote average over / base-pairs. Since sticky-ends and base- 
pairs near to both the ends of siRNA have large fluctuations compared to the middle part of siRNA, 
they are omitted for the calculation of 0. P(B) for siRNA while translocation inside various CNTs 
is plotted in Figure 6(c). We find that siRNA has to rotate in order to translocate inside CNT. 
Several mechanisms of DNA packaging propose the rotation of DNA during translocation [54- 
56]. By rotating, siRNA get reoriented with respect to the CNT inner surface which helps in 
overcoming local free energy barriers and helps in the translocation. Rotation is large for siRNA 
inside (20, 20) CNT at mM salt concentration and is least for siRNA inside (18, 18) CNT with 
P(B) having very strong peak at 132° and 35°, respectively. In the case of other CNTs, P(B) has 
more than one peak with comparable magnitudes indicating that those rotations are equally likely 
to be observed. Hence, the local free energy barriers are escaped by rotation of siRNA in order to 
translocate inside CNTs. 




(ID 
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IV. CONCLUSION 



To conclude, all atom molecular dynamics simulations with explicit solvent are used to inves- 
tigate the translocation and encapsulation of siRNA inside CNTs of various diameters, chiralities 
at various salt concentrations. After the translocation, the siRNA stays inside without any further 
movement. Free energy calculation using Umbrella sampling shows that siRNA gains in free en- 
ergy while translocating inside CNT and has to overcome large free energy barrier to exit from 
the interior of the nanotube. There is no translocation of siRNA inside (17, 17) CNT which has 
a diameter of 22.73 A. We find that the diameter of the (18, 18) CNT is the critical diameter for 
the translocation. Interestingly, siRNA translocates inside (20, 20) CNT but dsDNA cannot. This 
difference is due to (i) more interaction strength of uridine with CNT inner surface than that of 
thymidine [6, 7] and (ii) strong A-T base-pairing energy [6, 7, 47, 48] compared to A-U base- 
pairing energy. The translocation time x is decreasing with increasing diameter of the CNT with a 
critical diameter of 24 A. Inside CNTs of smaller diameters, vdW repulsion is larger which causes 
more structural deformations in siRNA. Free energy barrier for exit as well as translocation time x 
strongly depends on the chirality of the nanotube. The aromatic ring orientation in (25, 15) CNT 
also induces large structural deformations in siRNA and makes the siRNA and CNT interaction 
stronger than the (20, 20) CNT having same diameter. This makes the translocation in (25, 15) 
CNT faster compared to (20, 20) CNT. Among all the systems studied, siRNA has least defor- 
mations when translocated inside (20, 20) CNT. The stable siRNA-CNT hybrid can be used to 
deliver siRNA in RNAi applications. The kinetics and other thermodynamic analysis presented 
in this work allow us to get a microscopic understanding of the translocation process. By means 
of rotation, siRNA escapes from the local free energy barriers and translocate inside nanotubes. 
Modeling CNT with partial atomic charges on carbon atoms, functionalizing siRNA with polar 
nanoparticles for the efficient drug delivery systems are among exciting future perspectives in this 
emerging area. 
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(a)System after equilibration with water and counterions. 




(h)35 ns (i) 

FIG. 1 : (a) The initial system setup for simulation where siRNA-CNT hybrid was solvated with water and 
neutralizing Na + counterions. (b)-(i) Snapshots of siRNA inside (20, 20) CNT at 3 ns, 10 ns, 20 ns and 
35 ns during MD simulation in horizontal and vertical view with respect to CNT long axis h. Water and 
counterions were not shown in Figures (b)-(i) for clarity. The snapshots in this paper were rendered using 
VMD software [57]. 
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FIG. 2: The van der Waals interaction energy (j) (r) of siRNA that is arising due to the presence of CNT 
in the system. At t = ns, siRNA is outside CNT with two sticky-ends lying near to one end of the CNT. 
Hence at t = ns, the interaction of siRNA with CNT is very less. As time progress these sticky-ends would 
interact with CNT strongly as shown at t = 10, 15, 25 and 40 ns. 
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FIG. 3: (a) Binding free energy as a function of time with entropy and enthalpy shown in insets that are 
calculated using 2PT and MM-GBSA method, respectively, (b) vdW contribution to AG as a function of 
time. AG and (j) are increasing with time as siRNA translocates inside (20, 20) CNT and reaches constant 
after complete translocation. 
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FIG. 4: (a) Potential of the mean force F(£) of siRNA translocation inside (20, 20) CNT along the reaction 
coordinate £. (b) F{Q for siRNA translocation inside (20, 20) CNT at and 107 mM NaCl, (25, 15) CNT 
and dsDNA translocation inside (20, 20) CNT We define that the siRNA/dsDNA is translocated inside CNT 
when at least half of the base-pairs of siRNA/dsDNA are inside CNT from which it follows that £ < 59 A. 
The inset shows the minima of shifted F(%) which is inside CNT for siRNA and outside for dsDNA since it 
does not translocate. 
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FIG. 5: (a) The COM distance between siRNA and CNT (r com ) as shown in inset of supplementary Figure 
1(a) for siRNA inside CNT of various diameter and chirality. For the translocation, r com < L/2 implies 
r com < 43.5 A which is marked as horizontal dashed line, (b) Translocation time (x) of siRNA as a function 
of the diameter d of CNT using ff99 force field, x decreases linearly with increasing d. In the Figure, we 
have also shown x for siRNA in (20, 20) CNT using fflO force field as red square mark and siRNA in (20, 
20) CNT at 107 mM as diamond mark, (c) Number of close contacts (N c ) 
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FIG. 6: Structural deformation: (a) RMSd and (b) WC H-bonds of siRNA as a function of time. These 
structural parameters are fluctuating before complete translocation and reaches stable constant value after 
the complete translocation of siRNA. This also indicates that siRNA-CNT hybrid is stable after transloca- 
tion, (c) Probability of rotation angle, 6 of siRNA. 
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TABLE I: Summary of the simulation setup 



(m, n) CNT 


siRNA 


Box Dimensions (A 3 ) 


Na+ 


CI 


c {mM) 


WAT 


Total Atoms 
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(20, 20) 2880 


1396 


74 x 74 x 320 


44 





[NaCl] 


47619 


147177 


(20, 20) 2880 


1396 


74 x 74 x 320 


157 


113 


107 [NaCl] 


47395 


146731 


(19, 19) 2736 


1396 


73 x73 x318 


44 





[NaCl] 


45897 


141867 


(18, 18)2592 


1396 


72x72x319 


44 





[NaCl] 


44393 


137211 


(17, 17) 2448 


1396 


70 x 70 x 320 


44 





[NaCl] 


42690 


131958 


(25, 15) 2940 


1396 


74x74x318 


44 





[NaCl] 


47358 


146454 


(20, 20) 2880 


1398 


74x74x321 


42 





dsDNA [NaCl] 


47849 


147867 
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I. NON-EQUILIBRIUM VERSUS EQUILIBRIUM FREE ENERGY METHODS 



Transformation of binding energy from a function of simulation time t to a function of r com , 
center of mass distance between CNT and siRNA is carried out (See the schematic in Figure 1(a) 
for the pictorial definition of r com ). The transformation is done simply by monitoring r com (t). 
For the Umbrella sampling calculations, a new distance parameter called reaction coordinate 
is introduced. For simplicity in implementation without loss of generality and understanding, we 
calculate £, as the distance between the center of masses of 'far end' of CNT and first two Watson- 
Crick (WC) H-bonded base-pairs of siRNA near to the CNT (See schematic in Figure 1(a)). For 
the notations to follow, we use the term 'un-biased' for the MD simulation that are performed 
without any external biased forces and 'biased' for the MD simulations with external bias as the 
case of Umbrella sampling. In the un-biased MD simulations, the corresponding parameter for £, 
is ^un-biased- We plot ^un-biased as a function of r com (t) in Figure 1(a). Due to the effect of Umbrella 
biasing potential there could be minor structural deformations in siRNA leading to different values 
of ^ and ^un-biased for same value of r com in both the biased and un-biased simulations. Using the ^ 
values from Umbrella sampling, we plot ^ as a function of ^un-biased m the inset of Figure 1(b). By 
assuming (fitting) £, to depend linearly on ^un-biased; i-e., £, = «^ un _bi aS ed + c, we arrive at m = 1 .2022 
and c = -13.456. This calibration procedure helps us to calculate AG as a function of £ as shown in 
Figure 1(b) which is more relevant to compare with Umbrella sampling results. As can be seen in 
Figure 1(b), AG seems to have a global minima at ^ = 22 A where siRNA is inside (20, 20) CNT 
as the inset snapshot indicates. At ^ = 78 A or more, siRNA is away from CNT which is at peak 
of unstable configuration. Intermediate values of £, are also shown with snapshots. The limitation 
in calculating AG as a function of £, described above is that the system is not well sampled at each 
£, to have an ensemble averages since the translocation process is spontaneous. To improve the 
sampling at each £, to great level and obtain AG to great accuracy, we employ Umbrella sampling 
method as described in the previous section. Snapshots of siRNA in (20, 20) CNT are shown in 
Figure 2 at 40 ns in perspective display mode. 
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Figures 




(a) (b) 

FIG. 1: (a) Dependency of ^ U n-biased on r com from which £, can be calculated as ^ = 1.2^ un _bi aS ed — 13.456 
as shown in inset of Figure 1(b); where £ is the distance between the center of masses of 'far end' of 
CNT and first two Watson-Crick (WC) H-bonded base -pairs of siRNA near to the CNT as shown in inset 
cartoon, ^un-biased is the corresponding parameter in un-biased simulations, (b) AG as a function of reaction 
coordinate £, as denned in the Umbrella sampling for comparison. 



3 




(d) 

FIG. 2: Snapshots of siRNA after complete translocation inside (20, 20) CNT. These snapshots 
dered at 39 ns in perspective display mode. 



